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ABSTRACT 


The planetary obliquity plays a significant role in determining physical properties of planetary sur- 
faces and climate. As direct detection is constrained due to the present observation accuracy, kinetic 
theories are helpful to predict the evolution of the planetary obliquity. Here the coupling effect be- 
tween the eccentric Kozai-Lidov (EKL) effect and the equilibrium tide is extensively investigated, the 
planetary obliquity performs to follow two kinds of secular evolution paths, based on the conservation 
of total angular momentum. The equilibrium timescale of the planetary obliquity teq varies along with 
rz, which is defined as the initial timescale ratio of the tidal dissipation and secular perturbation. We 
numerically derive the linear relationship between teq and r; with the maximum likelihood method. 
The spin-axis orientation of S-type terrestrials orbiting M-dwarfs reverses over 90° when r; > 100, 
then enter the quasi-equilibrium state between 40° and 60°, while the maximum obliquity can reach 
130° when r, > 10*. Numerical simulations show that the maximum obliquity increases with the 
semi-major axis ratio a;/a2, but is not so sensitive to the eccentricity ea. The likelihood of obliquity 
flip for S-type terrestrials in general systems with a9 < 45 AU is closely related to m4. The observed 
potential oblique S-type planets HD 42936 b, GJ 86 Ab and 7 Boot Ab are explored to have a great 
possibility to be head-down over the secular evolution of spin. 


Keywords: planetary systems — planets and satellites: dynamical evolution — planet-star interactions 


1. INTRODUCTION 


'The planetary obliquity is simply defined as the an- 
gle between the planetary spin and normal orientation 
of the orbital plane. In our solar system, Uranus and 
Venus are well-known as retrograde planets with ex- 
tremely high planetary obliquity, where the obliquity of 
Uranus is 97?, whereas that of Venus is approximate to 
178? (Goldreich & Peale 1970; Dobrovolskis 1980). The 
primordial planetary obliquity provides a key clue to the 
understanding of evolution of our solar system, for ex- 
ample, the four inner terrestrial planets may have expe- 
rienced secular chaotic scenario of spin-axis orientation 
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(Laskar & Robutel 1993). For exoplanets, the evolution 
of planetary obliquity will directly affect the tempera- 
ture distribution on planetary surface, thereby leading 
to a variable seasonality of terrestrial planets (Gaidos & 
Williams 2004). Moreover, as to hot-Jupiters, high plan- 
etary obliquity can give rise to a severe atmosphere loss 
rate (Nikolov & Sainsbury-Martinez 2015). The obliq- 
uity can greatly influence the atmospheric escape rate 
of Earth-like planets orbiting late-type M-dwarfs (Dong 
et al. 2019), which plays a vital role in regulating surfi- 
cial habitability. 

As of today, 2M0122b is the only derived oblique ex- 
oplanet, which is believed to be a planetary mass com- 
panion around the 120 Myr host star at an orbital dis- 
tance of 50 AU. Thus, the diversity of planetary obliq- 
uity should be clearly understood due to planetary for- 
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mation. Bryan et al. (2020) measured three inclinations 
referring to the line-of sight, where to for the planet or- 
bit, ip for the planetary spin and i, for the stellar spin. 
The longitude of ascending node of the planet's equato- 
rial plane on its orbital plane is defined as Ay. According 
to the law of cosines, the most possible planet obliquity 
is estimated to be greater than 50?. So naturally arises 
a question: how does such high obliquity origin? 

Aside from indirect calculations through the spatial 
geometric configurations, a wide variety of observational 
methods were carried out to reveal the planetary obliq- 
uity. Kawahara (2016) analyzed the frequency modu- 
lation of the scattered light curve of a directly imaged 
exoplanet, which is induced by the planetary obliquity 
and orbital inclination. The modulation factor changes 
with difference attitudes of the spin rotations. Schwartz 
et al. (2016) further predicted that the planet's spin 
axis affects the time-resolved photometry, and east-west 
albedo contrast plays a vital role in constraining obliq- 
uity. Nikolov & Sainsbury-Martinez (2015) investigated 
the Rossiter-McLaughlin effect during secondary eclipse 
(RMse), and showed that for transiting exoplanets, the 
ratio of the ingress-to-egress radial velocity amplitudes 
are subjected to the planetary rotational rate and axial 
tilt. The synthetic near-infrared spectroscopic data were 
utilized to estimate the sky-projected spin axis orienta- 
tion and equatorial velocity of the planet. Currently, de- 
tection of planetary obliquity is limitedly known due to 
the spectroscopy accuracy. RMse is only effective when 
the host stars are brighter than K ~ 6 mag and large 
aperture telescopes (i.e., ~ 40 m) (Nikolov & Sainsbury- 
Martinez 2015). In this sense, kinetic theories are very 
helpful to understand the evolution of planetary obliq- 
uity in advance and its predictions will be further exam- 
ined by future high-precision measurements. 

Colombo (1966) proposed to describe the evolution of 
obliquity with the existence of a companion in the sys- 
tem, when the spin precession rate and orbital preces- 
sion rate become comparable, a resonance can occur to 
excite a larger obliquity. The spin axis motion is charac- 
terized as precessing around its orbital axis and the total 
angular momentum vector. The system configuration is 
specified by Euler angles (0, ¢, Y). The equilibrium solu- 
tions of the canonical equations of motion (Equation (6) 
in Fabrycky et al. (2007)) are defined as Cassini states 
(Colombo 1966; Peale 1969), where the planetary equa- 
torial plane maintains constant inclination to the plane 
of the ecliptic, and the spin axis remains coplanar with 
the normal to its orbit and the normal to the ecliptic 
plane. Total four Cassini states may exist for a given 
system, each with a preferred obliquity. 


Cassini states was adopted to explain moderately 
oblique exoplanet induced by planet-disk interaction 
(Su & Lai 2020) and the planet-planet interaction in 
multi-planetary systems (Su & Lai 2022a). The plane- 
tary obliquities can be maintained between 60° and 80° 
throughout the entire lifetime of Earth-like planets in 
the habitable zone of M-dwarf stars (Valente & Correia 
2022), which provides friendly temperate environment 
and conditions for habitability. 

In addition to the axial tilt of exoplanet, the spin-orbit 
alignment between planet and star plays a crucial part in 
the dynamical evolution of system. Storch & Lai (2014) 
intensively explored spin-orbit coupling effect between 
stellar spin and planet orbit, which leads to chaotic evo- 
lution of the stellar spin axis during Kozai cycles. Storch 
& Lai (2015) identified the origin of the chaos to be sec- 
ular spin-orbit resonances overlaps. Furthermore, Cam- 
pante et al. (2016) empirically constrained the angle be- 
tween a planet’s orbital axis and its parent star’s spin 
axis with the asteroseismology. Vervoort et al. (2022) 
indicated that the period and amplitude of the tilt of 
a planet’s rotational axis relates directly to long-term 
habitability of Earth-like planets. 

In this work, we primarily aim to investigate the plan- 
etary obliquity evolution of S-type terrestrials orbiting 
M-dwarfs in binary systems. The planets that simply 
revolve a single star in binary are referred to S-type or- 
bits, whereas those move around both stars are P-type 
orbits. As of May 2022, 154 binary star systems with 
217 planets were discovered, among which more than 
70 are S-type planets. For S-type systems, Kozai-Lidov 
model (Kozai 1962; Lidov 1962; von Zeipel 1910) is suit- 
able for a hierarchical three-body system with more in- 
tensive orbital precession when the perturbing object is 
more massive than the planet. In the secular evolution 
of a hierarchical triple body system, the perturbation of 
the third body exerts on a much longer timescale than its 
orbital period. The Kozai-Lidov mechanism for hierar- 
chical triple systems was studied under a wide variety of 
circumstances (Harrington 1968; Lee & Peale 2003; Naoz 
et al. 2013; Teyssandier et al. 2013; Tan et al. 2020). Li 
et al. (2014) explored Kozai-Lidov mechanism and char- 
acterized the parameter space that allowed large am- 
plitude oscillations in eccentricity and inclination under 
the test particle limit. For a system with an eccentric 
outer orbit, the octupole level terms in the perturba- 
tion Hamiltonian become significant, thereby triggering 
EKL effect. The planetary eccentricity increases approx- 
imately one, and the flips of orbital inclination can lead 
to a retrograde orbit (Naoz et al. 2013). Furthermore, 
the orbital flips of S-type planets were found to origin 
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from EKL (Huang & Ji 2022). Here we will explore 
similar flips of planetary spin-axis due to EKL. 

Tidal dissipation and the dynamical torque from the 
adjacent planets or disk both play critical roles in the 
spin evolution. Wang et al. (2019) adopted a compre- 
hensive scenario in the restricted three-body system, in- 
cluding the equilibrium tide and the third body pertur- 
bation, and showed that the obliquity librates for a long 
time or decays slowly down to zero. The planetary an- 
gular momentum, along with a positive dynamical loop 
consisting of orbital migration and obliquity tides, can 
well produce Ultra-Short-Period planets (Millholland & 
Spalding 2020). Formation of oblique super-Earths and 
several hot Jupiters can be further explained with tidal 
dissipation and spin-orbit capture (Su & Lai 2022b). 

In the real situation of a hierarchical system when the 
stellar companion is more massive than the perturbed 
body and the secondary's eccentricity is non-negligible, 
the non-restricted EKL mechanism will work to trig- 
ger the planet's evolution. Based on the conservation 
of angular momentum, in this work the mixed scenario 
between EKL and equilibrium tide is adopted to reveal 
a diverse secular evolution paths for planetary obliquity. 
The comprehensive evolution timescale of the obliquity 
depends on the masses, orbital elements of all bodies, 
and the tidal factors of the planet throughout the secu- 
lar evolution. The angular momentum transferred from 
the distant companion to the perturbed body can sus- 
tain or excite the planetary obliquity and contend with 
the tidal decay effect. 

To better understand the nature of planetary obliquity 
of exoplanets, we utilize N-body package MERCURY-T 
(Bolmont et al. 2015) to conduct more than 2000 nu- 
merical simulations that each run evolves at a timescale 
of roughly 107 yr to obtain two typical evolution paths 
of planetary obliquity corresponding two regions in the 
Hamiltonian level curves. Subsequently, we derive the 
relationship of equilibrium timescale versus the plane- 
tary mass m4, the initial orbit a, and the timescale ra- 
tio ri. In the parameter space of semi-major axis, we 
present the diagram of maximum obliquity to estimate 
the flip ratio of obliquity. Here we further investigate re- 
versal conditions of planetary spin axis involved in the 
EKL and tide effect simultaneously, which are employed 
to address the evolution of S-type planets, such as HD 
42936 Ab, GJ 86 Ab and 7 Boot Ab. 

This work is structured as follows. In Section 2, the 
theory of the mixed mechanism is described under the 
frame of central body equatorial coordinate, including 
theoretical analysis of the relationship between initial 
timescale ratio r; and the equilibrium timescale. Section 
3 describes evolution paths for general S-type terrestri- 


als, and provides the relationship between the maximum 
planetary obliquity, the equilibrium timescale and r;. 
The obliquity flip ratio is further calculated for a wide 
variety of initial conditions. Section 4 presents the secu- 
lar evolution of several potential oblique S-type planets, 
and predict the possible value of the equilibrium obliq- 
uity. In final, we summarize major outcomes. 


2. THE MIXED SCENARIOS 


In the secular evolution of triple body systems, the 
perturbation arising from the third body acts on a time 
scale much longer than the orbital period. If the pertur- 
bation of a circular outer orbit works, the Hamiltonian 
of this system can be expanded to the quadrupole level, 
which is referred to Kozai-Lidov mechanism as afore- 
mentioned. Here we adopt the Kozai-Lidov mechanism 
and the equilibrium tide to characterize the total evolu- 
tion for S-type planets. It is interesting that the rotation 
angular momentum, the planetary orbital angular mo- 
mentum and the stellar companion’s angular momentum 
are all participating to transfer to each other. The ge- 
ometry configuration of the system is plotted in Figure 
1, where the reference plane is the equatorial plane of 
the host star mo, mı and ma are the perturbed planet 
and the outer stellar companion. G'; and G2 are angu- 
lar momentum vectors of the inner and the outer orbits. 

is the total angular momentum vector of the orbital 
motion. N, is the angular momentum of planetary spin. 


2.1. Quadrupole Force vs Tide 


The classical Kozai-Lidov mechanism shows that the 
inclination and eccentricity of inner test particle oscil- 
late over secular evolution in a hierarchical system (von 
Zeipel 1910; Kozai 1962; Lidov 1962). With the secu- 
lar approximation, the inner and the outer orbits only 
exchange angular momentum, thus the semi-major axis 
of orbits do not change. When the SMA ratio a, = 
aı/a2 is a small parameter, the perturbation term of 
full Hamiltonian can be expanded as a power series in 
a, (Naoz 2016): 


k?momı k?mo (mo + mı) 


H= | 
2aı 202 


Ke. 2 rı j az j+1 
2 jm, (| Æ Er 
+ 2 ) aM; (2) (2) P;(cos ®), (1) 


where k? is the gravitational constant (with the mass 
unit of Mg and the length unit of au), rı is the dis- 
tance between mo and mj, ra is the distance between 
the center of mass of the inner binary and mz, Pj is the 
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Figure 1. Geometry of the hierarchical system. G, and à, are angular momentum vectors of the inner and the outer orbits. 


is the total angular momentum vector of the orbital motion. 


Legendre polynomial, ® is the angle between the vec- 
tors rı and ra, and the subscripts j = 1, 2 represent the 
inner and outer orbits, respectively. 

In the non-test particle approximation under the clas- 
sical Kozai-Lidov mechanism, the eccentricity and the 
inclination oscillate regularly in a well-defined timescale 
ti (Antognini 2015): 


16 a3 (1 ey? V mo + mi 


15 ar mk 


ly ~ (2) 
'This relationship was derived under the consideration of 
the equation of motion of w, argument of perihelion, by 
integrating between the maximum and minimum eccen- 
tricities. 

Meanwhile, we employ the constant time-lag equilib- 
rium tide at arbitrary rotation (Hut 1981; Leconte et al. 
2011; Wang et al. 2019) from the start of integration in 
the present work. Because the planetary eccentricity will 
be triggered to an extremely large value by Kozai-Lidov 
effect, which may result in unaccepable outcomes for 
the low-order expansion of eccentricity in the constant 
phase delay model. When the orbit plane does not coin- 
cide with the equatorial plane of the planet, Wang et al. 
(2019) used the Oxyz coordinate system to describe the 
equilibrium tide at arbitrary rotation and the equations 
of motion for the constant time-lag tidal model. 

We choose a central body coordinate system, and the 
stellar equatorial plane is treated as the reference plane. 
The origin is the host star, where x-y plane runs parallel 


p is the rotation angular momentum. 


to the equatorial plane of the star, z-axis is positively 
aligned with the rotation axis of the star. The Oxyz 
coordinate system is in line with the frame illustrated in 
Figure 1. When the tidal potential function expands to 
order 2 (Kozai 1965), the tidal acceleration is derived as 
(Wang et al. 2019): 


mo dT mi 


Fide = (Ff + Fo) + Fi), (3) 


momy 


5 . 
F, = mym [gme (Re) hy (1+ 3£r)] 


r? T 


Fo = 3878 (2) ka (9A -6) 7 


F, = 3@m (&)' ks), T[B cos(f +w) + C sin(f + w)], 
(4) 


A = cos Í cosi + sin I sini cos(Q — ©) 
B = cos I sin i — sin I cosicos(Q — ©) 
C = sin I sin(Q — ©) 
D =sinisin(Q — ©), 


(5) 


where w is the unit vector of orbital angular momen- 
tum, f is the vector from the star to the center of mass 
of the planet, and Ó is the unit vector along the orbital 
velocity. G is the gravitational constant, r is the dis- 
tance from the planet to the center of the star, ka is the 
2th order Love number, Å is the instantaneous orbital 
angular velocity. 7 denotes the time delay factor in the 
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Figure 2. The relationship between tr and t1 obtained from the theoretical equations. In panel (a)-(b), mı=0.2 Mọ, in (c)-(d), 
mi=1 Me. e» is set to be 0.2 in panels (a)-(c) and e2 —0.4 in (b)-(d). Dotted lines are the comprehensive evolution timescale 
of the planetary obliquity during the first output step of simulations. The comprehensive evolution timescale of the obliquity 
locally decreases when t1 > ti for ag > 10 AU and a; = 0.02 — 0.08 AU under EKL. 


constant time lag tidal model (Hut 1981). Q is the lon- 
gitude of ascending node of the planetary orbit, and i 
is the orbital inclination of the planet's orbit relative to 
the Oxyz system. Qp is the planetary rotational veloci- 
ties. I and ©, refer to the inclination and the longitude 
of ascending node of the planet's equatorial plane with 
respect to the Oxyz system, respectively. 

We can calculate equations of rotational motion based 
on the total angular momentum conservation: 


H + Hin + N = C, (6) 


where C is the total angular momentum vector, and the 


moment of inertia I, = mı R2r2, with r2 the square of 


radius of gyration (constant). Total angular momentum 


C and orbital angular momentum Hi, Ha can be calcu- 
lated by the initial orbital elements: ao, eo, io, Qo, Ip 
, Qp,o, Jo, Oo. In the evolution, the rotational velocity 
of planet Q, can be estimated with orbital parameters 
according to angular momentum conservation. Hence, I 
and O can be further derived through three components 
of the planetary spin vector. 

The evolutionary timescale of I is defined as (Wang 
et al. 2019): 


21,I,mıa 9 
t= m 2 a ( m 1)? 
3komg RŠ 7 
2mır?a$ 9 (7) 
= 2 ( e1)?, 
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where J, is the momentum inertia of the planet, T; is a 
typical tidal dissipation timescale, which keeps constant 
and can be calculated by: 


3 
, 
GmiTe 


T, (8) 
Tg is the time delay factor of Earth. In the following, we 
define the planetary obliquity as 0,, which is the angle 
between the orbital angular momentum vector and the 
rotational angular momentum vector (Figure 1). 

Note that Equation (1) and the secular perturbation 
timescale ty) in Equation (2) are defined in the Jacobian 
coordinate system (Naoz 2016). For a host star with 
mo — 0.1 Mo and the terrestrial planet with a mass of 
1 Ma, the center of mass of them will be inside the host 
star if the semi-major axis aı < 0.2 AU. The difference 
of coordinates between the Jacobian frame and the cen- 
tral body frame can be neglected. Here we employ the 
astrocentric coordinate as the system frame of reference. 

According to Colombo (1966), the spin precession rate 
and orbital precession rate is comparable to excite plan- 
etary obliquity. As the tidal dissipation timescale ty is 
much greater than txı induced by the Kozai-Lidov mech- 
anism. Blue and red solid lines in Figure 2 simply focus 
on the comparison of tjj and ty for a wide variety of 
initials m4, eo, a; and aa according to the theoretical 
equations, which is also considered to be a pair of initial 
conditions. Figure 2 indicates that when a; increases, 
the timescale ty can rise but tą) decreases. By contrast, 
ty grows when a» increases. For a, = 0.02 — 0.04 AU, 
the line of ty crosses fij at roughly 10% yr. The tidal 
dissipation timescale goes up along with the planetary 
mass, and tj is not so sensitive to the eccentricity e» 
of the secondary. In addition, Figure 2 exhibits that ty 
remains larger when ea = 0.2, ty and ty are equivalent 
(the intersections of two curves) at 5 x 10? — 3 x 10* 
yr. When ea = 0.4, the EKL timescale is smaller, and 
the intersection range of ty and tą is 4 x 10? — 3 x 104 
yr. For m4 = 1 Mo, the timescale of tide is larger in 
comparison with that of mı = 0.2 Ma, and the crossing 
point of tr and ty varies from 7 x 10? to 5 x 10% yr. 


2.2. Octupole Force vs T'ide 


Unlike the improved tidal triple body model (Wang 
et al. 2019), here we consider the combined effects of the 
octupole level perturbation and the equilibrium tide for 
non-restricted hierarchical systems. In such model, the 
total angular momentum is conserved, which consists 
of the orbital angular momentum of planet, that of the 
secondary star and the angular momentum of planet' 
spin. The three components of angular momentum are 
mutually coupling over secular evolution. However, for 


the restricted model, the total angular momentum of the 
planet is assumed to be constant. 

As Naoz (2016) pointed out, the EKL secular evolu- 
tion timescale is sensitive to ea of the secondary star, 
which is difficult to be directly calculated. Thus we 
numerically investigate secular evolution of the systems 
with N-body simulations. 

In Figure 2, red dots exhibit the comprehensive evo- 
lution timescale of planetary obliquity corresponding to 
different comparison between ty and tjj. Vertical co- 
ordinates of the red dots represent the mean evolution 
timescale of the planetary obliquity during the first in- 
tegration output step from t — 0 to t — ót. In the runs, 
the output step dt is set to be 100 yr. When the vertical 
coordinates of the dots fall down to be negative (gray 
lines) , the planetary obliquity goes down or evolves to 
be retrograde. For a relatively distant secondary star, 
0, will decrease initially with aa > 10 AU and a, = 0.02 
— 0.08 AU. 


3. SECULAR EVOLUTION OF THE OBLIQUITY 


Here the principal goal for this work is mainly to ex- 
plore the spin evolution of terrestrial planets, as the 
planetary obliquity plays a crucial role in habitability, 
therefore our study will cast light on spin evolution for 
habitable planets. Over two hundred S-type planetary 
systems were discovered so far, however, there are only 
10 binaries that own eccentricities, which those of the 
secondary stars are detected to be greater than 0.2. Note 
that the octupole force frequency is proportional to the 
eccentricity ea of the perturber under EKL scenario. 


3.1. Numerical setup 


'To extensively investigate the spin evolution for ter- 
restrial planets, here we adopt the standard equilib- 
rium tidal model and constant time lag model given 
by MERCURY-T (Bolmont et al. 2015) with numerical 
simulations. Here all runs refer to the N-body simula- 
tions. The constant time lag model consists on the as- 
sumption that the bodies under consideration are made 
of a weakly viscous fluid (Alexander 1973). MERCURY- 
T models the tidal forces between the star and the plan- 
ets but we neglect the tidal interaction between planets, 
then adds the forces in the acceleration of orbital mo- 
tion and the spin angular momentum. The tidal torque 
exerted by planet and the star is explicitly written in 
Bolmont et al. (2015). 

Here host stars are considered to be non-evolutionary 
in the simulations. Mass of the secondary star is set to 
be 0.2 Mo. The initials of aı are in the range 0.02 — 0.13 
AU, which keeps the orbital periastron of a rocky planet 
outside of the Roche limit of a host star (Liu et al. 2013). 
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For rocky planets, the lower value of Roche limit is 1.58 
Ro. The separation of the binaries az are given in [2, 5, 
10, 15, 20] AU. According to ty and txı in Equation (2) 
and (7), t1 can always be comparable with ty for each 
set of a2. To maintain the system stability, we choose 
eı = 0.2, e2 = | 0.2, 0.4] for terrestrial planets orbiting 
M-dwarfs, and set ea = [0.6, 0.8] for solar-type host star, 
respectively. The initial inclinations are given to be 24 
= 50° and ig = 0? with Qi = Q = 0°. 

Recently, the stellar spin was proved to has large im- 
plications on the dynamics of planetary systems (Becker 
et al. 2020; Chen et al. 2022; Faridani et al. 2023), 
our calculations demonstrate that the radius of M-dwarf 
with age of 700 Myr and a mass of 0.1 Mo will decrease 
from 0.1238 Ro to 0.1236 Ro. The terrestrial planet's 
orbital angular momentum variation due to the stellar 
spin evolution over 2 x 107 yr is less than 4 x 106, 
and the final planetary obliquity differs from that of the 
non-evolving model by ~ 1°. Thus we adopt the non- 
evolving M-dwarf host model, which will not bring about 
significant influence to our results. 

In the non-evolving host body model, the time delay 
factor in the constant time lag tidal model is set as Te 
— 698 s (Bolmont et al. 2015). For terrestrials, we select 
the love number kz = 0.305, the radius of gyration r? = 
0.3308, and the planetary radius Rp = Rg. The initial 
planet rotation period is set to be 24 hr, and the initial 
value of 0, is assumed to be 10° in the simulations. 

The former study revealed that general relativity 
(GR) may have influence on triples’ evolution, which 
may result in resonant-like dynamics (Naoz et al. 2013; 
Liu et al. 2015) or destabilize the existing resonance 
(Hansen & Naoz 2020). We conduct additional runs 
by considering GR effect with the same initials as Fig- 
ure 4. For a terrestrial planet at a; = 0.04 AU with the 
GR and non-GR models, the difference in final rotation 
period is 0.5 hours and the final planetary orbital an- 
gular momentum differs by 6 x 107°. To speed up the 
integration, here we simply adopt the non-GR model in 
the simulations. 

In this work, we conduct over 2000 simulations that 
each evolves for 10 Myr (~ 10% ta) with MERCURY- 
T. When all runs reach 107 yr, the simulation is ceased 
once the equilibrium state of planetary obliquity occurs 
after 10" yr, which is the stopping condition. 


3.2. Typical Evolution Paths 


According to the origin of Cassini States described in 
this work and Equation (8) in Su & Lai (2022b), we first 
plot the Hamiltonian level curves in the phase space of 
cos 0 — ó to theoretically describe the spin evolution fea- 
ture, @ is the angle to represent the precessional phase 


of spin vector about the orbital angular momentum vec- 
tor. We find the dynamical structures similar to Su & 
Lai (2022b) as shown in Figure 3, including up to three 
kinds of typical oscillation modes, which refer to three 
kinds of Cassini States. For region I in Figure 3, the 
obliquity suffers excitation below 90?. 0, experiences a 
larger oscillation amplitude and the spin-axis flip with 
0, > 90° in region II. After flip, 0, stably librates be- 
tween 40? and 60?. In region III, the planetary obliquity 
is initially stirred up to 0, > 90°, finally evolves to the 
stable retrograde spin orientation. 

The Hamiltonian level curves indicate that the planet 
spin evolution dynamics is greatly affected by the ini- 
tial conditions, the equilibrium obliquity corresponding 
to Cassini State II that decreases with the precession 
ratio of the orbit and spin-axis. Based on the theoret- 
ical results from the Hamiltonian level curves, we then 
mainly investigate the typical oscillation cases of region 
I and II, to summarize the relationship between initial 
conditions and the final evolution feature. 

Figure 4 presents a typical set of obliquity evolu- 
tion with respect to diverse initial orbits for terrestrials, 
where a, of the planet varies from 0.02 to 0.2 AU. Panel 
(a)-(d) follow the evolution path of region I in Figure 3, 
and panel (e)-(f) follow the region II. We also perform 
the time-series evolution of inclination and eccentricity 
in Figure 5 to show the effect of the EKL strength, and 
the excitation effect of EKL is greater when the planet 
and the secondary star are closer. The eccentricity can 
be excited to 0.7 and the orbital inclination decreases to 
lower than 40° when a; = 0.1 AU. The oscillation ampli- 
tude of eı and i, are increasing with a4. The variation 
of semi-major axis is too little during the evolution to 
capture the orbit decay, thus we could not see the no- 
ticeable circulation. 

The evolution results can be further distinguished by 
rz, as shown in Figure 6. Here we perform the numerical 
results with respect to four groups of m; = 0.2 Ma, 1.0 
Me and ea = 0.2, 0.4. With the increase of r+, the 
maximum obliquity also increases. Figure 6 shows that 
the planetary obliquity without flip when r; « 1, while 
the obliquity of Earth-mass planet could flip when r; 
ranges from 1 to 10°, and 6, will always flip when r, 
spans from 10? to 10°. 

It is noteworthy that when 10 < r < 10°, Omar goes 
up to a local bulge and then declines to a valley. It is 
speculated that the extreme planetary eccentricity origi- 
nates from EKL mechanism when r; gets larger, and the 
tidal effect is dramatically enhanced owing to a much 
closer pericenter to the host star. Under the combined 
scenarios, the planetary obliquity undergoes an inhib- 
ited growth. On the other hand, we further compare 
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Figure 3. Three examples of Hamiltonian level curves for the spin evolution of terrestrial planets around the M-dwarf. The 
planetary semi-major axis ranges from 0.05 AU to 0.1 AU, and the precession frequency ratio noted as 7 by (Su & Lai 2020) is 
selected as [0.12, 0.45, 1.2] respectively. The different Cassini states are marked out with labels I, II, III. 
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Figure 4. Evolution of the planetary obliquity or S-type terrestrials with a mass of 1 Ma and semi-major axis smaller than 
0.2 AU, other initial conditions are set to be mo = 0.1 Mo, mı = 1 Ma, ma = 0.1 Mo, ai = [0.02, 0.04, 0.045, 0.05, 0.1, 0.12, 


0.14, 0.2] AU, aa = 5 AU, eı = 0.4, e2 = 0.2, iı = 50 deg, i? = 0 deg. The obliquity can be triggered to retrograde and the 
time to reach the obliquity equilibrium state is much sensitive to the planetary initial position. 
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the obliquity evolution due to es. When e» rises to 0.4, 
the more obvious obliquity excitation occurs, and the 
maximum 6, enters the flip region with smaller logr. 
For r; > 10%, the maximum obliquity converges to ~ 
130°. 

In total numerical simulations, we calculate the evolu- 
tion paths for terrestrial planets with mı = 0.2 ~ 8 Ms, 
where the general evolution trend is consistent. When 
mı grows to 1 Mg , the local maximum obliquity of 
the bulge ascends to 100° with aa = 20 AU. Additional 
simulation results reveal that when m increases to 4 
Ma and 8 Ms, the local maximum obliquity can arrive 
at 105? and 110?, respectively. It can be inferred that 
more runs may locate at the region of flip with larger 
planetary masses for terrestrials around M-dwarfs. 

We carry out additional N-body simulations with 
MERCURY-T to study more general S-type systems 
with an extremely higher ea. Currently, there are over 
five S-type binary systems that are observed to have an 
eccentric outer companion with ea > 0.5, thus we can 
further investigate the evolution of 0, in these cases. 
The simulations show that the terrestrials orbiting M- 
dwarfs can be scattered or captured by the secondary 
with ea = 0.8. Furthermore, we conduct a couple of runs 
to explore terrestrials orbiting solar-type stars in com- 
parison with those about M-dwarfs, where 55% of them 
can survive under the comprehensive effect of EKL and 
tidal dissipation. 

Aside from typical oscillation cases presented in Fig- 
ure 4, we also find chaotic evolution behaviour. Storch & 
Lai (2014) demonstrated the stellar spin-orbit angle can 
have a chaotic region, and similar chaotic behaviour of 
the planetary obliquity exist in our work. These chaotic 
cases always occur when ag = 2 AU and es > 0.6, the 
extremely excitation of eccentricity leads to the instabil- 
ity of the system. Comparing with the stellar spin-orbit 
angle chaotic behaviour, the chaotic results of planetary 
obliquity in our work origin from spin-orbit coupling ef- 
fect between the planet spin and the planet orbit, rather 
than between the stellar spin and the planet orbit. 


3.3. The Equilibrium Timescale 


In Section 3.2, 0, remains a quasi-equilibrium state for 
over 10” yr prior to a decrease to zero, and the timescale 
of entering the equilibrium obliquity is denoted as teq- 
We then calculate the ratio between the orbital period 
and rotation period P. over the entire evolution. For all 
terrestrial cases with mı = 1 Ma and a, = 0.1 ~ 0.2 
AU, P, will fall down to 6 at the time of teq, which seems 
to be in a high-order spin-orbit resonant-like state. The 
total angular momentum of the planet, the summation 
of Gi and Ny, also drops down to a constant. We infer 


that the planetary angular momentum is trapped in the 
high-order spin-orbit resonant-like state. 

Su & Lai (2022a) showed that the equilibrium obliq- 
uity ranges from —90? to 90? for Cassini State I, II and 
IV, whereas the obliquity varies from 160? to 180? for 
Cassini State III. The related precession ratio of spin 
vector and orbital angular momentum vector is in (107? 
—10!). Here the equilibrium 6, performs the oscillation 
state consisting of Cassini State II- evolves between 40? 
and 60?, which agrees with the typically evolution in 
our simulations. Winn & Holman (2005) addressed that 
Cassini state II is the most favorable configuration to 
maintain a significant obliquity. In this study, we adopt 
a non-restricted model, the initial planetary mass is set 
to be 0.2 - 8 Ma, then we perform in a larger parameter 
space of r, spanning from 107? to 107. Our findings un- 
veil a new equilibrium value of the maximum obliquity 
Omax close to 130° when r; > 10* (Figure 6). 

In the equilibrium stage, the initial planetary orbit re- 
siding in the habitable zone of M-dwarf is not shrunk to 
approach the Roche-limit, thus the spin-axis orientation 
remains relatively stable, thereby helping maintaining 
the planetary habitability. 

Figure 7 shows the contours of equilibrium timescale 
teq in the parameter space aı-aa, where teq evenly dis- 
tributes in the dimension of a, but varies with az and 
e2. Next, we observe that teq is always less than 10° yr 
for a; = 0.09 AU (see the dashed lines in each panel). 
For a? < 20 AU, teq is approximately proportional to 
the initial orbit of planet, which primarily governs tidal 
timescale. For aı = 0.02 AU and a; > 0.1 AU, the obliq- 
uity equilibrium timescale keeps unchanged for the vari- 
ations az at same a1. When the semi-major axis meets 
0.04 < a, « 0.08 AU, aa = 20 AU, the tidal timescale is 
equivalent to the secular perturbation timescale. 

We further investigate the relationship between teq 
and the timescale ratio r; as shown in Figure 8, where 
Figure 8 (a)-(b) shows the planetary mass of m; — 0.2 
Mag while mı = 1Mg in Figure 8 (c)-(d), with respect 
to eg = 0.2 and 0.4, respectively, for each panel. We 
then utilize a maximum likelihood (ML) function and 
numerically optimize the relation. The fitting model is 
logteg = Mm x (logr;) + bart. The fitting results of teq 
and tj/tyj are provided in the panel legend. When az 
= 2 or 5 AU, the tidal torque and the EKL effect are 
intensive for (tr/t3) < 10?, the equilibrium timescale 
of the obliquity much arises from the tidal dissipation, 
thus the blue dots deviate to the lower end. 


3.4. The Obliquity Flip Ratio Raip 


In Section 3.1, we describe two typical evolution paths 
of obliquity, and cases follow region II can always trig- 
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Figure 5. Evolution of the eccentricity and inclination for S-type terrestrials with a mass of 1 Mẹ and semi-major axis smaller 
than 0.1 AU, other initials same as Figure 4. The eccentricity can be excited to 0.7 and the inclination decreases to below 40° 
when a; = 0.1 AU. The oscillation amplitude of e1 and i1 are increasing with a1. When aı > 0.1 AU, the upper and lower limit 


eı and i, in the evolution are not sensitive to a1 any more. 


ger the flip of spin-axis. We extend the distance of the 
second star to 45 AU in this section, and investigate the 
obliquity flip ratio in the parameter space of aı-aa. In 
Figure 9, over 700 simulations are conducted for S-type 
terrestrials with mı = 0.2 ~ 8 Ma, accompanied with 
a constant perturber of 0.2 solar mass, in which each 
panel is plotted with a grid resolution of 10 x 18. The 
flip ratio is calculated by the fraction of the simulations 
that undergoes flip. 

The flip ratio in Figure 9(a)-(b) are 0.49, 0.55, 0.67, 
0.72, respectively. We conclude that the reverse of the 
spin axis can be very common for S-type terrestrials 
with a distance companion. When the planetary mass 
increases, the flip ratio increases with the same initial 
semi-major axis, thus the flips get easier with higher 
m/m». This will provide the clues to predict what kind 
of exoplanets can evolve to be head-down, like Venus in 
our solar system. We also change the eccentricity of the 
companion, but the flip region of obliquity is not much 
sensitive to ea. For the Earth-mass planet, the flip ratio 
spans from 0.52 to 0.56 when e2= 0.2 — 0.6. 


Table 1. Parameters for binary stars in 
three potential oblique S-type systems 


HD 42936 GJ 86 r Boot 
aAB(AU) 1.22 20 45 
CAB 0.594 0.4 0.91 
ma(Mo) 0.87 0.8 14 
mp (Mo) 0.08 0.49 0.4 


4. APPLICATIONS TO POTENTIAL OBLIQUE 
S-TYPE SYSTEMS 


In this section, several S-type planets are reported to 
figure out whether they rotate with the oblique attitude 
under the combined scenarios of EKL and equilibrium 
tide. Considering the timescale of secular evolution, we 
select the terrestrial planet HD 42936 b and the Jupiter- 
like GJ 86 Ab and 7 Boot Ab as case study. In Table 
1 and 2, we summarize essential parameters for three 
potential oblique S-type systems. 
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Figure 6. Distribution of the maximum obliquity Omas in the dimension of timescale ratio r+. The semi-major axis a2 ranges 
from 2 to 20 AU. Regions with or without flip are also filled with pink and green respectively, with 90? as horizontal boundaries. 


Table 2. Physical parameters for planets in three 
potential oblique S-type systems. 


HD42936b GJ 86 Ab r Boot Ab 


ap(AU) 

P(day) 
ep 

ip (deg) 

wp (deg) 

Mp (Mz) 


0.066 
6.67 
0.140 


80.6 
0.008 


0.110 
15.77 
0.046 
266.0 
4.01 


0.046 
3.31 
0.079 
45 
218.4 
4.13 


4.1. HD 42936 b 


HD 42936 (DMPP-3) is the closest binary harboring 
a S-type super-Earth HD 42936 b (Barnes et al. 2020). 
The primary star HD 42936 is a solar-like star with mass 
of 0.87 Mo at a distance of 48.9+ 0.6 pc, while the semi- 
major axis of the secondary star is simply 1.22 + 0.02 
AU, and the eccentricity is relatively high with eg = 
0.594. Such a compact and hierarchical system naturally 
brings secular perturbation scenario and tidal evolution 
into our mind. Barnes et al. (2020) concluded that the 
EKL mechanism, as one of the most likely formation 
models, can drive HD 42936 b move to current orbital 
configuration, in which the eccentricity of planet can be 
initially stirred up, and then shrink to a circular orbit 
when the tidal timescale is comparable to or shorter than 
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Figure 7. The equilibrium timescale map in the parameter spaces of semi-major axis. In panel (a)-(b), the planetary mass 
mı = 0.2Me, e» = 0.2 and 0.4, in panel (c)-(d), the planetary mass changes to mı = 0.2Ma, and ea is same as panel (c)-(d). 


tij. Hence, HD 42936 is a good example to explore the 
spin evolution with our mixed mechanism and reveal the 
dependence of the spin-axis evolution tendency on the 
initial conditions. 

HD 42936 b was discovered by radial velocity, prob- 
ably moving at an inclined orbit. Thus in this work, 
the initial semi-major axis of this planet is assumed to 
be 0.06 — 0.12 AU, and the initial inclination ranges 
from 10° to 170°. When 40? < ip < 140°, HD 42936 
b will collide with the primary star or be scattered out 
of the system, before that the obliquity 0, can be ex- 
cited to 150°. For the survivals, 0, can be triggered to 
the maximum value of 80?. Figure 10 shows the obliq- 
uity evolution for HD 42936 b with the initial prograde 
orbits, while Figure 11 stands for those results for the 
retrograde orbits at beginning. 


4.2. GJ 86 Ab 


Aside from terrestrials in binary systems, we again 
investigate the spin evolution of hot-Jupiters with the 
mass ranging from 0.3 to 8 My. The simulation results 
show that the obliquity can be triggered to retrograde, 
and the maximum obliquity can arrive at 160°, being 
indicative of that the equilibrium obliquity may be sen- 
sitive to the planetary mass. 

GJ 86 is a nearby S-type system at a distance of 10.9 
pc from the solar system, and hosts a giant planet GJ 
86 Ab (Queloz et al. 2000) at the orbital distance of 
0.11 AU, the outer companion is a white dwarf in the 
eccentric orbit with ag = 21 AU, ea = 0.4 (Mugrauer 
& Neuhäuser 2005; Lagrange et al. 2006). The orbital 
elements and masses of both companions are measured 
(Zeng et al. 2022) with the observations from radial ve- 
locity and high angular resolution imaging along with 
the absolute astrometry data from Hipparcos (Perry- 
man et al. 1997) and Gaia (Gaia Collaboration et al. 
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Figure 8. The relationship between teq and the timescale ratio ri. Panel (a): mi = 0.2Mo, e2 = 0.2; Panel (b): mi = 0.2Me, 
e2 = 0.4. Panel (c): mi = 1Me, e2 = 0.2; Panel (d): mi = 1Me, e2 = 0.4. In each plot, the initial semi-major axis of the 
secondary star ranges from 2 AU to 20 AU. Dots plotted with different markers are numerical simulations results. The solid 
lines are the estimated relationships fitted by the maximum likelihood method. 


2016; Brandt 2018). The current distance between the 
binary is 21 AU, whereas the mutual closest approach 
may reach 9 AU (Zeng et al. 2022), being indicative of 
that the gas-giant may origin from the disk truncation. 
The secular perturbation also participates in the obliq- 
uity evolution due to close separation of the binary. Our 
simulations suggest that GJ 86 Ab have a large possibil- 
ity to follow obliquity evolution path II, and maintain 
the head-down rotation attitude. 

GJ 86 is selected as one of the candidates for CHES 
(Closeby Habitable Exoplanet Survey) mission (Ji et al. 
2022), which will observe nearby solar-type stars to hunt 
for terrestrial planets in the habitable zones at an ultra- 
high resolution via astrometry. The future observations 


will extensively provide the orbital elements and true 
mass of GJ 86 Ab, which helps place constraints on the 
dynamics of spin evolution. 


4.3. r Boot Ab 


The giant planet orbiting r Boot is also one of the 
best-known exoplanets in the nearby stars from our so- 
lar system. The planet 7 Boot Ab was first discovered 
in 1996 (Butler et al. 1997). Brogi et al. (2012) uti- 
lized radial velocity measurements to determine an or- 
bital inclination of i, = 44.5 + 1.5 ° and a true plan- 
etary mass of 5.95 + 0.28 My. They presented the at- 
mospheric characterization of this non-transiting planet 
and showed the planetary temperature distribution de- 
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Figure 9. The spin-axis flip diagram. Each panel has specific initial parameters mı and ea, (a): mi = 0.2Mg, e» = 0.2, (b): 
mi = 0.2Me, e» = 0.4, (c): mı = 1Mg, ea = 0.2, (d: mi = 1Me, e» = 0.4. The colour bar represents the range of the 
maximum value of 05 over 100 Myr, whereas those in blue denote cases without spin-axis flip. The flipping regions are filled 


with yellow-green and orange. 


creases towards higher altitudes. Table 1 and 2 summa- 
rize the essential parameters for 7 Boot. For a, = 0.046 
AU, aap = 45 AU, we still believe the perturbation from 
the outer companion can excite the orbital eccentricity 
and spin obliquity of r Boot Ab, as the eccentricity ea 
of the secondary remains extremely high. Our numeri- 
cal results indicate that 0, follows the evolution path II, 
and the maximum value can reach above 100°, then falls 
to the equilibrium stage before 100 Myr (Figure 12). 


5. CONCLUSIONS AND DISCUSSION 


In this work, we explore the spin evolution of terres- 
trial planets as the planetary obliquity plays a vital role 
in habitability. As the binary systems of stars that har- 
bor planets are very common, thus secular perturbation 
from the secondary may greatly lead to a diverse plan- 


etary evolution. Here we adopt a comprehensive model 
composed of the octupole level secular perturbation and 
equilibrium tide to study dynamical evolution of S-type 
planets’ obliquity, which may provide essential clues for 
elucidating spin evolution for habitable planets. 

We first calculate the comprehensive evolution 
timescale of 0p, which relates very much to the Kozai- 
Lidov timescale tjj and tidal timescale tj. We further 
plot the Hamiltonian level curves of the spin dynamics 
to theoretically derive the evolution paths of terrestrials' 
planetary obliquity, and also conduct numerical simula- 
tions to extensively explore the evolution of planetary 
obliquity induced by the coupling effect due to EKL res- 
onance and the equilibrium tide, the planetary obliquity 
equilibrium state varies as a function of r;. In addition, 
the maximum 6, during evolution can be distinguished 
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Figure 10. The obliquity evolution for HD 42936 b of the 
orbital inclination 4; —10?, 20°, 30°. 
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Figure 11. The obliquity evolution for HD 42936 b of the 
orbital inclination 4; —140?, 150°, 160°. 


by rı. Cases in region I presents no extreme excitation 
of the obliquity, and 0, will stably oscillate between 40? 
and 60? in region II, the planetary obliquity has excita- 
tion with flip. When es ascends to 0.4, more obliquity 
excitation emerges, and 0, with five sets of aa enter the 
flip region earlier. For r; > 10*, the maximum planetary 
obliquity evolves to nearly 130°. 

According to the numerical results in the parameter 
space of semi-major axis a, and ag, we present dia- 
grams of the equilibrium timescale and maximum obliq- 
uities. The equilibrium timescale teq evenly distributes 
in the dimension of a,, and the timescale is always 
smaller than 10° yr when a4 < 0.09 AU. For az < 20 
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Figure 12. The obliquity evolution for r Boot Ab of the 
orbital inclination i; = 45°. 


AU, teq is approximately proportional to the initial or- 
bit of the planet, which mainly determines the tidal 
timescale. Then we derive the linear relationship be- 
tween teq and the timescale ratio rs with the maximum 
likelihood method, logteg is approximately proportional 
to logr; and the slope of fitted lines rises with m. 

Moreover, we estimate the ratio of obliquity flip in 
the extensive parameter space ag « 45 AU. We con- 
clude that the reverse of spin axis appear to be common 
for S-type terrestrial planets accompanied with a distant 
massive companion. When the planetary mass ascends, 
the flip ratio grows with the same initial semi-major axis, 
thus the flips become easier for a larger mı /ma. To em- 
ploy our comprehensive model in the observed systems 
and predict the evolution of potential oblique planets, 
then we further numerically investigate three cases of 
HD 42936 b, GJ 86 Ab and 7 Boot Ab, to understand 
their evolution paths of 0p. 

Aside from the planetary spin-orbit flip, the flip of 
the star-spin-orbit angle induced by the EKL mechanism 
were noted in the literature previously (Naoz et al. 2011, 
2012; Petrovich 2015a; Lai et al. 2018; Stephan et al. 
2018), and the original high stellar spin-orbit obliquity 
will be damped for the cooler star (Lai 2012; Dawson 
2014). As above-mentioned, Storch & Lai (2014) inten- 
sively studied the spin-orbit coupling effect between the 
stellar spin and the planet orbit, which gives rise to the 
chaotic evolution of stellar spin axis during Kozai cycles. 

Effects of tidal dissipation and stellar spin-down will 
also influence the final distribution of spin-orbit mis- 
alignment angles of hot-Jupiters. Storch & Lai (2015) 
explored the origin of this chaotic stellar spin behaviour, 
they identified secular spin-orbit resonances and the res- 
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onance overlaps are responsible for the chaos. Key pa- 
rameters including the adiabaticity parameter, the ratio 
of the Kozai-Lidov nodal precession rate and the stellar 
spin precession rate are proved to dominate the degree 
of chaos. In our research, similar chaotic behavior of the 
planetary obliquity also exist in several simulation cases. 
It was also demonstrated that before being destroyed 
by their stars either Roche-limit crossing or engulfment 
during stellar expansion, there could be diverse results 
on the final stellar-spin orbit angle (Petrovich 2015a,b; 
Stephan et al. 2018; Angelo et al. 2022). 

It was recently highlighted that chaotic tide plays a 
key role in the evolution of planetary systems (Wu 2018; 
Vick et al. 2019). In this case, in eccentric orbits, the 
planet's orbital energy is converted into internal fluid 
energy via tidal stretching and compression at perias- 
tron passage, causing the orbit to circularize over time. 
For Earth-like planets, a multilayered internal structure 
would increase the efficiency of tidal dissipation and af- 
fect the climate and habitability of the planet. For 
Jupiter-mass planets, increasing the rate of planetary 
orbital energy loss accelerates the tidal dissipation pro- 
cess and faster entry into the proposed equilibrium state. 

With the improvement of observational accuracy, the 
space-borne telescopes, such as JWST (James Webb 
Space Telescope) (Gardner et al. 2006), TESS (Tran- 
siting Exoplanet Survey Satellite) (Ricker et al. 2015), 
PLATO (Planetary Transits and Oscillations of stars) 


(Rauer et al. 2014), ARIEL (Atmospheric remote- 
sensing infrared exoplanet large survey) (Tinetti et al. 
2021) and CHES (Closeby Habitable Exoplanet Survey) 
(Ji et al. 2022), will release a bulk of detailed informa- 
tion about the planetary three dimension orbital con- 
figuration, the planetary mass, the presence of atmo- 
sphere, and spin orientation, thereby further identifying 
the stable oblique planet in habitable zone. There are 
prospects for constraining exoplanet obliquities in the 
coming years, such as using high-resolution spectroscopy 
to obtain the projected rotational radial-velocity of the 
planet (Snellen et al. 2014; Bryan et al. 2018) based on 
the Rossiter-Mclaughlin effect during the second transit 
eclipsing. The planetary obliquity is also an important 
indicator for the seasonality of habitable terrestrials and 
the atmosphere loss rate of hot-Jupiters. This can be 
further examined by upcoming observations. In future 
work, we will further explore the origin of these chaotic 
behaviors, the stellar spin evolution, and the connection 
between planetary obliquity and stellar obliquity as well. 
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